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We consider the 3-pt function (i.e. the bispectrum or non-Gaussianity) for stochastic backgrounds 
of gravitational waves. We estimate the amplitude of this signal for the primordial inflationary 
. background, gravitational waves generated during preheating, and for gravitational waves produced 

' by self-ordering scalar fields following a global phase transition. To assess detectability, we describe 

I how to extract the 3-pt signal from an idealized interferometric experiment and compute the signal 

to noise ratio as a function of integration time. The 3-pt signal for the stochastic gravitational wave 
background generated by inflation is unsurprisingly tiny. For gravitational radiation generated by 
purely causal, classical mechanisms we find that, no matter how non-linear the process is, the 3-pt 
correlations produced vanish in direct detection experiments. On the other hand, we show that in 
QQ ' scenarios where the B-mode of the CMB is sourced by gravitational waves generated by a global 

phase transition, a strong 3-pt signal among the polarization modes could also be produced. This 
may provide another method of distinguishing inflationary B-modes. To carry out this computation, 
we have developed a diagrammatic approach to the calculation of stochastic gravitational waves 
' sourced by scalar fluids, which has applications beyond the present scenario. 
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I. INTRODUCTION 



Great strides are being made in both the technological and theoretical aspects of the direct detection of astrophysical 
^ ' gravitational waves. Ground based interferometers like LIGO 1] have achieved their design sensitivities and are in the 
process of being upgraded to even greater precision 0]- Design of a space-based detector, LISA is well underway 
with a pathfinder mission scheduled for launch in mid-2011 and a possible launch date for the full mission in the next 
decade or so. Further, planning for the next generation space-based gravitational wave detectors, BBO and DECIGO, 
has begun 0, Q . These detectors are specifically designed to search for cosmological stochastic gravitational waves 
(SGW). Beyond direct detection experiments, large scale B-mode polarization of the cosmic microwave background 
[a| can be sourced by a spectrum of long wavelength gravitational waves at last scattering. The Planck mission will 
constrain a tensor to scalar ratio on the order of r '--^ 0.1 while work is well underway on proposals ,7i] that will probe 
to r~ 0.01. 

A stochastic background of gravitational waves can be generated in a variety of ways. Unresolved point sources such 



as neutron star or black hole binary systems generate a stochastic background in the confusion limit [8]. Quantum 
, fiuctuations of the metric during infiation [1, [lO, [HI [13] are amplified on super Hubble scales, generating a stochastic 
background whose amplitude directly probes the energy scale of infiation. Following inflation, explosive par ticle 
production associated with a phase of pre/reheating would also produce a stochastic background p^. [T^. [Tsl. Il6l ITtI. 
^ ■ [l8, 19, 20, 21, 22]. Further, SGW can be generated by phase transitions 23, 24, 25, 26, 27, 28|, bubble collisions 



I [2 3 . |30j or even more exotic processes involving warped extra dimensions 31 1 . 

Recently, models have been proposed where a scale-invariant spectrum of SGW, mimicking the inflationary spec- 
trum on subhorizon scales, is generated by the relaxation of a disordered scalar matter^ field with a white noise 
spectrum on sup er horizon scales. Such a white noise spectrum can be laid down by a global phase transition event, 
for example [26l . [27l [3^ . More generally, non-equilibrium processes in cosmology produce scale-dependent spectra of 
gravitational waves. This should not be a surprise: violent motion of large masses generates a non-trivial quadrupole 
moment. 

In the study of cosmological density perturbations in large scale structure and in the Cosmic Microwave Background, 
such higher correlations of perturbations are called non-Gaussianities [33| - any Gaussian field is completely described 
by its power spectrum, or 2-pt function.^ Indeed, while inflation is expected to produce a highly Gaussian spectrum 



^ In this paper we will deal exclusively with scalar matter although this is not a necessary condition. 

^ For a completely Gaussian field, all correlation functions of an odd number of fields vanish and all even correlations can be written as 
products of the power spectrum. Such truncated correlation functions are usually called disconnected, see for example j34| . 
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of SGW akin to that of the spectrum of density perturbations [35|, cosmological SGW foregrounds from active sources 
are expected to be highly non-Gaussian. However, these processes occur on characteristic scales which are smaller 
than the size of the horizon. We present a simple argument that any causal process, i.e. one which operates inside 
a single post-inflationary horizon volume, which generates gravitational waves predicts a vanishing 3-pt correlation 
in our detectors. Such processes include preheating (which we will discuss in detail as an example below), bubble 
collisions, and all of the other processes mentioned above. Inflation on the other hand predicts correlations among all 
modes at all scales. 

In addition to the present work, the study of non-Gaussian features in SGW spectra has received little attention: 
Drasco and Flanagan investigated the statistics of popcorn noise 13611 . while Seto has suggested the use of the 4-pt 
correlator to study intermittent bursts [13, Ull . Racine and Cutler [39| investigated deviations from Gaussianity of 
unresolved galactic white dwarfs binaries and suggested that it is small due to the large number of sources. 

This paper is organized as follows. In Section |lT] we discuss sources of 3-pt correlation functions and their prospects 
for detection - in Section III Al we consider the 3-pt signal generated by quantum effects during inflation while in 
Section III Bl we estimate the expected amplitude of the 3-pt function from active scalar sources. We consider two 
scenarios, global phase transitions and preheating, and present some preliminary analytical results. The details of 
the calculations are left to the Appendices. In Appendix |^ we construct the 3-pt estimator of SGW, given a set of 
detector data streams. In Appendix |B] present a diagrammatic method of calculating general iV-pt SGW correlation 
functions from linear scalar sources which will have applications beyond the present scenario. Finally, we conclude in 
Section Uni 



II. SOURCES AND OBSERVATIONS 



Cosmological SGW can be divided into those sourced by initial quantum fluctuations (namely those laid down 
during inflation) and those sourced classically by non-zero quadrupole moments. Gravitational waves from inflation, 
at the moment of creation, behave like free fields with Gaussian initial conditions and hence possess a vanishing 3-pt 
correlation. Non-zero 3-pt correlations are formed when these fields interact gravitationally (see the left graph in 
Fig. [T]) and carry information about the gravitational coupling strength H/Mp. Inflationary gravitational waves are 
correlated on all scales due to their creation during an epoch of accelerated expansion. These correlations are laid 
down as the scales leave the horizon and are frozen until reheating. On the other hand, actively sourced gravitational 
waves are generated by gravitational hremsstrahlung, and carry information about the quadrupole moment of the 
source (see the right graph in Fig. [1]). These gravitational waves are laid down as each scale enters the horizon, and 
are only correlated on scales comparable to the horizon size at the time they were created. We consider two speciflc 
models, preheating and the global phase transition scenario of Jones-Smith et al. [26l. [33|. 

A 3-pt correlation of gravitational waves is produced during inflation via the right graph in Fig. [TJ However, this is 
highly suppressed as each internal scalar line is a copy of the scalar power spectrum, supressing the graph by a factor 
oiO{V{kf) - 10-30. 



A. Inflation 



Even if inflation is exactly de Sitter, the spectrum of SGW it generates will be non-Gaussian, due to gravitational 
self-interaction. At leading order, this non- Gaussianity is sourced by a 3-pt interaction term, represented diagram- 



matically in Fig. [T] (left graph). The amplitude of this process is calculated using the "in-in" formalism^ 41 1. 

Each external GW leg is 0{H/Mp) due to canonical normalization of the graviton, while the coupling term is 
purely gravitational, and also 0{H / Mp). Thus, after freeze out, gravitational waves outside the horizon will possess 
a non-zero 3-pt correlation function, with amplitude (H/Mp)'^. 

In any model of inflation in purely Einstein gravity, the 3-pt correlation function of gravitational waves at freeze-out 
is [H 

(/i^(k)/i^'(k')/.^"(k"))inf = (-K + + ^'J + + ^] (2nrS{k + + k") 



X 



— ^(e^l,(k)ef,.,(k')eil, 0^%,^,,,,), (1) 



%{kk'k"f ^ " ^ ' " 



^ In general, these amplitudes will also be corrected by both scalar and GW loops, but they will be subdominant - see for example [40ll . 
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FIG. 1: Non-trivial 3-pt correlations of gravitational waves can be generated directly from graviton-graviton interactions (left 
graph) or indirectly via interactions with scalars (right graph). Both processes are always present, however, direct interaction 
dominates the 3-pt correlation during inflation while for scalar sourced gravitational waves this process is highly suppressed 
relative to the loop. 



where i/* is the Hubble rate at horizon crossing, K — k + k' + k" , and the tensorial structure is = k'^Sji + k"^5ii + 
k^5ij. In this work, since the exact details of the polarization sum in Eqn. ([1]) do not concern us, we replace it with 
the appropriate powers of k from dimensional analysis. 

To calculate the present amplitude of the inflationary 3-pt correlation (Eqn. ([T|)) we evolve the metric perturbations 
from horizon re-entry using the gravitational wave transfer function.^ For modes reentering the horizon during 

— 1/2 

radiation domination, A^^u = Zeq (Hq/2h f), replacing k = 2nf^l. We find 

•^(/, ,/= /")today = ^(/, /, /")inf X J Jl Jll "3/2 ' ^"^^ 

where T{f^ /', /") is defined by 

(/i(/, ti)h{f', n')h{f", n")) = T{f, /', f")S{fn + fn' + f'n"). (3) 

Here, z^q « 4 x 10"^ is the redshift of matter-radiation equality and Hq is the Hubble constant today. As detailed in 
Appendix \^ the total integration time T of a direct detection experiment can be chopped into "chunks" of identical 
duration AT ^ /*, where /, is roughly the frequency of minimum noise. We can then construct the 3-pt estimator per 
chunk and, since the noise is uncorrelated across chunks, the total signal to noise scales as VM, where M ~ T/AT 
is the total number of chunks: 

SNR = .hape=< (1=1)* (g)' (^^g.) ^ ( ^,(;.)^,|(^|(,.K^;)3 ) >< W 

where A/ is the width of the window in frequency space. Here we collect the frequency integral into a dimensionless 
quantity we call "shape" ^ 

shapes f'^ dcos0;, f^'''^^ fdff''df'f"^df" (^) ' L I"):J-^5{f - Vf'^ + f"^ + 2/'/" cos 0^), (5) 

where 9^ is the angle between n' and ft" . The shape encodes where the 3-pt signal has support in the 3 dimensional 
parameter space (/, /', /") - different gravitational wave generation mechanisms will, in general, produce different 
shape functions. The shape can be integrated numerically, but we note that the filters pick up triangles that are 
roughly equilateral, and approximate the shape simply as 

1088 , if^ + Ap) A/3 

shape = ——TT- ^— (6) 

9 (/2 _ Apf U 



* In reality the GW transfer function is much more complicated l42l . however, the corrections are 0(1) and so for this work we neglect 
them. Also, one might worry that the different reentry time for each mode will induce a phase shift in the 3-pt correlation today. 
However, it is easy to show that for modes reentering in the radiation era, the phase shift in the signal 3-pt is an overall constant and 
hence can be set to zero. 

^ This nomenclature follows standard convention in the studies of CMB non-Gaussianitios, where the the 3-pt is similarly integrated over 
all possible triangles to produce a single measure, see for example O . 
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Given a total integration time T, we obtain M = T/(10 x f~^) independent observations, including a fudge factor of 
10 to reduce edge effects from chopping the signal. The total integration time expected for a 90% confidence detection 
then scales like (assuming identical, time independent and constant detector noise) 

T^fl'N^if,)iAf)-^. (7) 

In other words, as in the detection of the SGW power spectrum, low noise at a low target frequency /* is advantageous. 

An instrument like BBO is designed to detect the inflationary power spectrum, but it is highly unlikely that we 
will be able to detect the 3-pt correlation - this is a consequence of the fact that inflation is highly Gaussian. For a 
typical BBO/DECIGO detector with N{f^) ~ lO-^^Hz-^/^ at - 0.1 Hz and A/ ~ 0.033 Hz, and assuming GUT 
scale inflation i7inf/Mpi ^ 10~^, T ^ 0.27 x 10^h~^ years where Hq = lOOft. km/s/Mpc. It may appear that, since T 
scales with /* (assuming /* ~ A/), if we were to increase the sensitivity of BBO by one order of magnitude 

or move /* down a decade in frequency, we would be able to detect the 3-pt correlation on a reasonable timescale. 
However, this result scales with {Hini / Mp\)~^ . Unless the energy scale of inflation is near the upper limit currently 
allowed by WMAP [35| measurements of the cosmic microwave background, the inflationary 3-pt function will be well 
out of reach of direct detection experiments. 



B. Active Scalar Sources 



Beside the stochastic gravitational waves generated by the amplification of quantum fluctuations during inflation, 
gravitational waves can be generated classically via violent motion of mass-energy. For example, a period of turbu- 
lent cosmological evolution, say during preheating, will result in the generation of large gradients and the copious 
production of SGW. 

"Active sources" refers to sources which are physically moving. We focus on the results here and present a general 
review of scalar sourced gravitational waves in Appendix |B] We will see that for all causally generated, actively 
sourced, gravitational waves, the 3-pt function vanishes in direct detection experiments. Gravitational waves from 
preheating will not be detectable in the CMB, since the long wavelength modes necessary to generate the temperature 
and polarization anisotropics are not produced, but global phase transition models continuously sources horizon scale 
gravitational waves. These modes will affect the CMB anisotropics on scales smaller than the horizon (in real space) 
at last scattering. In particular, we will show that the source of the polarization modes in the CMB for this model 
has large 3-pt correlations (i.e. it is highly non-Gaussian). This means that the polarization bispectrum of the CMB, 
i.e. both the (BBB) and the (EBB) statistics will be non-trivial for this model, and could provide a powerful test 
for non-inflationary sources of the polarization modes. 



1. Preheating and other causal mechanisms 

During inflation, quantum fluctuations of the metric are amplified on super Hubble scales by the accelerated 
expansion. Following inflation, in most models, the inflaton decays and reheats the universe. The first stage of this 
process, preheating, is dominated by an explosive and non-perturbative production of highly inhomogeneous, non- 
thermal fluctuations of the inflaton and the other fields coupled to it. The inhomogeneous decay of the inflaton and 
the turbulent phase that follows it are inevitably accompanied by the production of gravitational waves. This topic 
was flrst discussed by Khlebnikov and Tkachev [l^l and went largely untouched for almost ten years. ^ Recently a 
flurry of papers has appeared with more accurate numerical simulations improving and expanding on earlier work 
[13, EH E, di lJ9j 20, 21j, and exploring the scaling relationship for this signal. 

The stochastic gravitational wave spectrum frorn preheating has a peak amplitude at a physical scale that depends 
only on the energy density at the end of inflation [ij] , 

^^Klfd- (8) 

This can be expressed as a frequency via 



/.6xW^Hz, (9) 



See also [3. 



5 



where Hf. is the Hubble rate at the end of inflation. To obtain a peak at / = 0.1 Hz requires He/Mp\ ~ 10~^^, or 
inflation ending near the TeV scale. The peak amplitude of the preh eating gravitational waves is (largely) independent 
of the frequency, and is estimated to be ^gwh^ ^ 10^^^ today p^ . 

Since preheating is a completely causal process, only modes that are within the horizon are excited. If kp is the 
physical scale of preheating then 

> 1. (10) 

Negligible amounts of gravitational waves are generated outside the horizon, so we do not expect super-horizon 
correlations. Modes within each Hubble patch are expected to be highly correlated at the 3-pt level, we expect 
gravitational waves to be uncorrelated across patches. Since the signal in our detectors is the sum of contributions 
from a large number of uncorrelated patches, a simple application of the central limit theorem implies that their 3-pt 
correlation function will be highly suppressed relative to the 2-pt function. Specifically, there will be approximately 
TV ~ (He/Ho)'^ ~ 10*° patches for TeV scale inflation (and even more for GUT scale inflation). The 3-pt correlation 
function for preheating must be suppressed by (He/Ho), i.e. the square root of the number of different causal 
disconnected patches in the sky, relative to the 2-pt function.^ We thus expect that SGW backgrounds from preheating 
would be observed to be highly Gaussian today. If such a spectrum is detected by direct detection experiments there is 
the intriguing possibility that this highly gaussian source of SGW may be used as a backlight to probe the foreground 
structures. We postpone a discussion of this possibility to a future publication. 

This argument for Gaussianity does not rely on any of the details of preheating and thus extends to other SGW 
that are sourced on subhorizon scales by active processes in the early universe. 



2. Global Phase Transitions and the CMB 



It has been suggested that a global phase transition in the early universe can actively source a scale invariant 
spectrum of stochastic gravitational waves on large scales, mimicking that of inflation [26j. It has also been suggested 
that such a process may also mean that the detection of large scale B-mode polarization of the CMB would not 
be a unique "smoking gun" signal of inflation Using causality arguments, Baumann and Zaldarriaga [i^l have 
shown the polarization signal from an actively sourced spectrum of gravitational waves would be distinct from that of 
inflation. The usual B-niodes are defined non-locally in terms of the Stokes parameters Q and U, and are not obliged to 
vanish for scales outside the horizon. However, Baumann and Zaldarriaga construct a real space correlation function 
which respects causality. For inflation, this correlation function has features on scales larger than the size of the 
causal horizon at last scattering that cannot be present for any mechanism which causally generates gravitational 
waves during the standard, post-inflationary period of the universe. We will see that, even if a such a scenario can 
completely mimic inflation at the level of the 2-pt function, they are vastly different at the level of the 3-pt function. 
This may provide an additional test for determining the inflationary origin of the CMB B-mode polarization signal. 

We consider the model of [1^ (see also [l^] , the formalism of which we adopt for this paper) and present only the 
key details here. The reader is referred to the original papers for additional information. A field, $(x, t), in the vector 
representation of 0{N) is in an initially symmetric state with zero vacuum expectation value (vev), (^(x, i)) = 0. 
The field is maintained in this state by a quadratic potential which possibly arises from thermal corrections, or from 
a coupling to the inflaton. As the universe evolves, thermal corrections become negligible or inflation ends via a 
tachyonic instability which causes the potential to evolve into a mexican hat. The field obtains a non-zero vev, 
(^(x, i)) — V, by "rolling" to the true vacuum state. On scales larger than the horizon size, the direction in field 
space of the vector $(x,i) is uncorrelated and there is gradient energy associated with the — 1 Goldstone modes, 
p ~ (9$)^. As these modes enter the horizon, the scalar field aligns itself and some of the gradient energy is radiated 
into gravitational waves. 

On large scales (or low energies) the field is confined to the vacuum manifold, 0a (^i ^) — '^'^i ^^'^ the dynamics 
of the A^ — 1 Goldstone modes are well described by the non-linear sigma model. Furthermore, in the large A^ limit. 



This same argument can be applied to compact binary sources; while individually each source is highly non-Gaussian, the central limit 
theorem tells us that as long as each individual source is uncorrelated the sum will be Gaussian. In principle, although preheating 
patches "know" about the inflaton potential and hence is correlated in some way, the fact that the process is likely to be highly chaotic 
means that this knowledge is rapidly lost and hence we do not expect the patches to be correlated. Nevertheless, this is not a given - 
recently [45ll argued that for some classes of SUSY-inspired models of inflation, large non-Gaussian spikes of the curvature perturbation 
can be laid down over superhorizon scales which preserve this memory, perturbations which may source a super-horizon spectrum of 
SGW. We are extremely grateful to Richard Easther and Lam Hui for pointing this out to us. 
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the model is soluble Initial conditions (at r = r,) are assumed to be white noise on super horizon scales with 
vanishing power on subhorizon scales, corresponding to the field being initially aligned on these scales 

iUKrM^\r.)) ^ I 6.^r3(2.)3i#.(k + k') , k,. << 1 
[ , k-q^ > 1. 

It is easy to see why a strong correlation between the 2-pt and 3-pt correlation functions is expected from such a 
source: both are generated from the same interaction term (Fig[3]) at 1-loop. From dimensional analysis, M^h ~ {4>4')j 
and in Fourier space 

k'ihl) - m)f , k^hl) ^ {{cj^Df. (12) 
Via Wick's theorem for a linear source (f>, it follows that 



1. (13) 



In other words, we should expect that the 2-pt and 3-pt (and indeed, any higher point) functions are equally impor- 
tant.® Since gravitational waves source both E and B modes, we expect that such a mechanism sources non-trivial 
{BBB) and {EBB) correlations on the CMB sky. As both such correlations are vanishingly small in the standard in- 
flationary scenario, measurement of these correlations will be a smoking gun for a non-standard source of polarization^. 
We leave the details of the construction of the 3-pt polarization correlation functions to future work. 

Let us now return to the detailed calculation, one can skip right ahead to Eqn. (|21[) if one is not interested in the 
technical details. For power law expansion in conformal time r, a oc (/3 = 1 for radiation domination, 2 for matter 
domination), one finds (0a(k, r)(/)b(k', r')) — (27r)3(5at5(k -t- k')F(fc, t, t') where 

F(k r r') v' r(/3 + l/2)r(2/3 + 3/2) 3,, Ji+,(fcr) J.^kr') 

The 3-pt function induced by this source in the co- located detector approximation is then given by Eqn. (|B17p . 
which together with the above "propagator" (technically a stochastic average) for the scalar field gives 



(/i,,(k,/)/.,,(k',/')/^H(k",/")) (15) 
„„(k)0,fe,op(k')OH,,r(k")^(k + k' + k")^^^iV (^67r2^: 



r/Ti dTn dTo n.( )t? sin ffc-Ti InYxT Wi? sin ffc'T-olnYxo Wi? sinffc"xol / 

(27r)3 

J2(|k' + p|t3) Jadk' + p|t2) J2{pTi) J2{PT2) J2(|k - p|ri) J2(|k - pirs) 



X J dTidT2dT3a{Ti)T^am[kTi]a{T2)T^sin[k'T2]a{T3)T^sm[k"T3] J J^^^™-^' 
^{poPpip- k)q{p- k) 

+ {p- k)o{p- k)pPqPr 



(|k' + p|r3)2 (|k' + p|r2)2 (pri)2 {pr2Y (|k-p|Ti)2 (|k-p|T3)2 
■h{K + pIts) J2{\k" + p|t2) J2(|k - Pin) Jadk - p|r2) J2(pri) 72(^X3) 

(|k"+p|r3)2 (|k"+p|T2)2 (|k~p|Ti)2 (|k-p|r2)2 {pT,Y ipT,Y 



Here C'ij_fc;(k') is the transverse traceless projector, defined in Appendix [BJ An exact analytic evaluation of this 
expression is difficult. In principle there is no obstacle to numerical integration, but for our purposes it is sufficient 
to approximate it as follows: 

• Work in the equilateral limit, |k| = |k"| — |k"|, $1 -f Jl' -f ft" = 0. Most of the power is produced at horizon 
crossing, so we expect that the signal will be strongly peaked on equilateral shapes. 

• Work in the long wavelength limit kr < 1, kr' < 1 and fcr" < 1 for all times between t* and Tend — 1/fc so that 
sin(fcr) « kr etc 

• We neglect the angular dependence of |k" + p| and |k — p| 



In the language of the scalar bispecturm, complete correlation between (ff) and (CCC) would mean ~ 10^. 

In fact, this process is present even during inflation - gravitons are sourced by bremsstrahlung of the inflaton field itself, but it is easy 
to show that the contribution to the inflationary 3-pt from this process is highly suppressed. 
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• We can then use asymptotic expansions of the Bessel functions. In the range p < min(l/Ti, 1/t2, I/T3) we use the 
small argument expansion of the Bessel function. In the region min(l/Ti, I/T2, 1/T3) < p < max(l/Ti, I/T2, 1/T3) 
we distinguish between large and small argument expansions and finally in the range max(l/ri, l/r2, 1/T3) < 
p < 1/t, we can use the large argument expansion of the Bessel function. 

There are two limits in which we can expand the Bessel functions, large and small argument expansions, these are 
Under these approximations, we can write 

{h,J{■k,f)h,k{k',f')hk^{k"J")) (17) 

« -40,^.„„(k)0^.fe,„p(k')0fe.,,.(k'')^(k + k' + k'')^^^iV (^evr^^^MElZZ^^ /™„(fc). 
The integral is 

/"■"^(A:) = 3!2i/3r!3/2 f'^'^ dn H dT2 H dr^ k^T^T^Tl ( /"'^"' dp/i"''' + /'^"'dp/2""' 

J-n, J-q, Jri, ^Jo Ji/ti 

^^jrnnopqr^ / dp I^'^'P'^'X (18) 

1/t2 J I/T3 ^ 

where the factor of 3! accounts for permutations of the integration variables, {ri, T2, T3}, and we have used the scale 
factor during radiation domination a « Hq^/U^t, consistent with ao ~ 1 today. The integrals are 

^ (2^ / - ^)'(^^ - (19) 



(27r)326 TT (pTl)5 TT (pT2)-'5 

_ 1 2 cos^ (pn + f )^ 2 cos^ (pr2 + 2 cos^ (pr3 + f ) / ™ n o p. _ _ 
- (27r)3,r (pri)5 tt (pr^^ tt {pr.f J dU p p p p {p k) [p k) 

To perform these integrals, we replace the cos^(pt) terms by their value averaged over a few cycles, 1/2. The angular 
integrals can be done in the usual way. In the equilateral limit 

Oy,m„(k)0,fe,op(k')(^fc^,9r-(k") y"dr!p™pVpW = ^, 

fc°FO,,-„„(k)0,fe,op(k')OH,gr(k") Id^r^rtPf = (20) 



20 1 / w V fHoV^^ ' 



Finally on subhorizon scales we obtain 

(kkk)..(k + k' + k", ^ [^^) *(k + k' + k"). (21) 

The subscript A denotes the equilateral limit, iliad^^ — 4.15 x 10~^ is the radiation density today and N is the 
number of components in the scalar field, which is taken to be ~ 4. 

We can calculate the 2-pt function in the analogous way, using the propagator above in Eqn. (jB18|l the integral is 

{h.,,{k,T)h,,{k',T)) = Ar(^6^2^lW^|£|I/^^ 0,,,„„(k)(i|^^''dTi^''dr2a(Ti)sin(fcri)a(r2)sin(fcr2) 

J2{pri) J2{PT2) J2(|P - k|ri) J2(|p - k|T2) 

7 (27r)3 (pTiY {pT2y (|p-k|ri)2 (|p-k|r2)2 
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In the same approximation as above, we find for modes inside the horizon 



Nk^ VMpJ k 



(/..,(k,r)/^,,(k',r))~-^( — ) (^^V^) 5(k + k'). (23) 



Then, writing 

k^{hhh)A=CNL{k^{hh)f^^, (24) 
where Cnl is a dimensionless constant, for this theory we have 

20 



7r3/2\/]V 



0{1). (25) 



for N — 4. For inflation Cnl ^ 0{H/Mp\), while for this model the 3-pt function is as important as the 2-pt function. 
Note that Cnl is independent of the scale at which the process is occurring. This means that, even if a global phase 
transition can completely mimic inflation at the level of the power spectrum, it produces a 3-pt which is distinct from 
and much larger than that of inflation. The existence of such a large 3-pt correlation function relative to the 2-pt 
correlation function suggests that if a B-mode signal consistent with r ^ 0.01 was observed, by constructing the 3-pt 
estimator we might determine whether its origins were consistent with inflation. 

We also point out that this 3-pt would vanish in direct detection experiments for precisely the same reasons as in 
the preheating case. At first it may seem that since the gravitational waves arc continuously sourccd that we might 
be able to evade the simple argument above based on the central limit theorem. However, the power at each scale, 
k, is sourced as that particular mode enters the horizon. Since direct detection experiments are sensitive to scales 
on the order of the size of the solar system, the gravitational waves detected in these experiments will be primarily 
composed of radiation which was emitted when the horizon scale was on the order of the size of our solar system. 
This means that these gravitational waves will again look almost completely Gaussian by the central limit theorem. 



III. CONCLUSIONS AND FUTURE OUTLOOK 



In this paper we considered the properties of the 3-point statistics of cosmological gravitational waves from both 
inflationary and non-inflationary "active" scalar sources. For the latter, we consider gravitational waves from pre- 
heating driven turbulence at the end of inflation and from self-ordering scalar flelds following a global phase transition 
in the early universe. 

Introducing a "3-pt correlation parameter" Cnl, we write 

k^{hhh)^=CNL{k^{hh)fl'^, (26) 

where the a here denotes the equilateral limit of the bispectrum. During inflation, metric fluctuations (gravitons) 
start in a purely Gaussian state and interactions with other fluctuations are highly suppressed by the amplitude of 
the observed scalar spectrum. This means that the leading order effective three graviton interaction is the tree level 
interaction, and any departure from Gaussianity is highly suppressed. For GUT scale inflation, Cnl ~ ^finf/-^pi ^ 
10-<^, and Cnl measures the the scale of inflation directly, which provides a consistency check when compared to the 
SGW power spectrum. Unfortimately, and unsurprisingly, we find that even a BBO/DECIGO class detector cannot 
detect an inflationary 3-pt signal. 

On the other hand, for SGW from scalar sources, Cnl measures the strength of the effective gravitational wave 
interaction, i.e. it is mediated by scalars. For these sources there is no reason to expect the higher point functions to be 
suppressed relative to the 2-pt function. In highly inhomogenous phenomenon that follows any period of cosmological 
turbulence, we estimate that Cjvl ~ C(l) ^ Hinf/Mp. 

We investigated two different scalar sources of SGW. The first model is that of gravitational waves produced in the 
era of preheating after inflation. During preheating, gravitational waves are produced with a characteristic wavelength 
kp ~ He where Hg is the Hubble scale at the end of inflation. While the gravitational waves are highly correlated at 
this scale, since '> Htoday, such sources appear as uncorrelated patches in the sky. That is, these patches appear 
to our detectors as uncorrelated point sources. The gravitational radiation we observe in our detectors is the sum of 
gravitation waves coming from all directions, and is thus composed of gravitational waves from many patches. By 
constructing the correlator of three time streams, we are integrating the estimator {h{n, f)h{n', f')h{n", /")) over the 
entire sky. Clearly the 3-pt signal from such a source will be highly Gaussian via the central limit theorem - ironically 
even more so than the inflationary signal. Nevertheless, the gravitational radiation from processes such as TeV scale 
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preheating is expected to be easily seen by future direct detection experiments. The high level of Gaussianity of these 
stochastic backgrounds then presents the intriguing possibility of using them as a probe of the intervening matter 
distribution. 

The second model we consider is gravitational raditation from self ordering scalar fields following a global phase 
transition proposed in this context by Jones-Smith, Krauss and Mathur 26]. In this model, a scale invariant spectrum 
of SGW is generated by self ordering scalar fields at the horizon. We show that the correlation parameter Cnl ~ 
0(1) in this model and a large 3-pt correlation is expected. Nevertheless, since these gravitational waves are being 
continuously sourced at the horizon, the central limit argument above still applies and the 3-pt signal today will be 
small to vanishing. 

One potential place where a non-Gaussian gravitational wave signal may still be detected is in the polarization 
measurements of the CMB. The gravitational waves sourced by the global phase transition mechanism are completely 
correlated, which means that we expect that as well as the usual B-mode spectrum, {BB), they will also give rise 
to non-trivial (BBB) and (EBB) correlations, and could provide an independent discriminant of non- inflationary 
sources of B mode polarizations. 
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APPENDIX A: THE 3-PT ESTIMATOR FOR DIRECT DETECTION EXPERIMENTS 

In this Appendix, we calculate the correlation of three streams of data from a direct detection experiment such 
as LISA. A detector output is a scalar stream of data S{t) as a function of time, where the data measure the direct 
strain of the detector. The data stream consists of a signal s{t) and an inherent noise n{t) 

S{t)^s{t)+n{t). (Al) 

In general we are deep in the noise dominated regime, \n\ 3> |s|, (s) = (n) = 0, and for two detectors in different 
locations, the noise is assumed to be localized, (riiri2) <SC {n-in-i), (71,2712). In the measurement of a stochastic power 
spectrum, we can correlate two detectors to extract the signal out of the noise. In this case, the signal to noise ratio 
for the stochastic power spectrum is 



SNR^-^i^, (A2) 



where (S'12) is the expectation value of the signal 



5i2= / dt^l dt2Si{h)S2{h)Q{h,t2), (A3) 

J-T/2 J-T/2 



and Q{ti,t2) is a 2-pt filter function. The denominator of Eqn. (|A2p is the root mean square noise \/ (-^12): where 
A^i2 — S12 — (S'12). Since the noise is locahzed, (S'12) — J dti J dt2(si(^i)s2(^2))W^(ii, ^2)- 

By choosing a filter with a moving window, Q{ti,t2) = Q{ti — t2), we see that the signal to noise ratio in Eqn. 
(|A2[) scales like Vt, where T is the total integration time. The signal in the numerator is correlated and hence scales 
like T while the denominator is uncorrelated, and so increases like a one dimensional random walk, Vt. 

With three detectors, we can construct the analogous estimator for the 3-pt using a filter W{ti,t2,t3), 

T/2 

dhSi{h)S2{t2)S^{h)W{tiMM)- (A4) 

T/2 





.T/2 


.T/2 


S123 — 


/ dti 


/ dt2 / 




J-T/2 J 


'-T/2 J- 
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Localization of the detector noise means that 

(^123)=/ dt^l dt2 I dh{si{h)s2{t2)si{h))W{ti,t2M)- (A5) 



.T/2 


.T/2 .J 


/ dh 


/ dt2 


/-T/2 J 


-T/2 J- 



The noise here is N123, = S123 — {S123), and the signal to noise ratio is SNR ~ {S123) / y/ (-^123) • 

To observe a signal at a frequency /* we must make a measurement that is at least AT ^ 1//* in length. Then, 
instead of making one long measurement over the total time T, we chop our signal into "chunks" of length AT ^ 1//*, 
where /* corresponds to the characteristic minimal noise frequency of a given detector, i.e. Na{f*) = minimal, where 
the subscript a here labels the detector. We then have T/AT chunks of identical measurement of the signal (si23)m, 
where M labels the M-th chunk. On the other hand, the noise is uncorrelated across chunks and, as noted above, 
increases like a one-dimensional random walk. The signal to noise per chunk is then SNRa/ ~ {S123) m / \J -^1231 where 
(<S'i23)m is the expectation of Eqn. (|A4p after an integration time AT ~ 1//*- 
The total SNR is then given by 

(si23)Af 



SNR = ^ 'PJI X VA/ cx VT. (A6) 

That is, the signal to noise ratio of a 3-pt correlator also scales like VT. This means that a detection is inevitable as 
long as we integrate for long enough (and the 3-pt is actually non-zero). 

We now need to relate the estimator in Eqn. (|A4[) above to the predicted 3-pt signal, (x, ti)hjk{x' , i2)/ifci(x", t^)) 
where the Roman indices run over the 3 spatial dimensions. We can expand any massless, transverse-traceless tensor 
mode, hij (x, t) , as 



hij 



(x,i) ^JdnJ /2d/^/^(0,/)e-2-/(*-"--)e^(r2), (A7) 



where k = 27r/r2 , k = 27r/ and the speed of light c — I. Our convention differs from that of the gravitational wave 
community, who usually absorb the into the amplitude (see for example [isil . Iso!]). but conforms with standard 
analytical techniques used to compute higher order correlation functions. The polarization tensors are normalized via 

ef^{k)e^'{k)^26^^', (A8) 

where polarization is indexed with A. The explicit form of the polarization tensors can be found in (for example) [50| . 

A gravitational wave is a spin-2 tensor field propagating through space at the speed of light. To detect a gravitational 
wave, it has to interact with a detector such as an interferometer or resonant mass. Each detector, labeled a, has a 
pattern tensor D^J , which measures its response to a passing GW as a function of its geometry. We define the detector 
pattern function as 

Dl^et^{n) = F^{n). (A9) 

A detector located at x^ then produces a data stream S'a(i,Xa) — Sa(t,x.a) + na{t,x.a) in response to a passing GW, 
where Sa{t, x^) is 

Sait,^a)= j fdfd(lY,h^{f,a)e-^^^^f^*-^---^F^{n). (AlO) 

With three data streams 5'i(ii), 6*2(^2) and 5*3(^3), and imposing the filter function W{t,t' ,t"), the expectation of 
the estimator (5'i23) in Eqn. (jA5P is related to the GW 3-pt function by 



(5i23) - / dt I dt' [ dt"dndn'dn"fdfrdff"^df" V {h^{f,n)h^' {f',n')h^" if",n')) 

J J~T/2 J -T/2 . 



A.A'.A'' 



xF^{tl)F^\tl')Ff {tl'')e-^^'^f^'''^-^'^+f''^''-'^'-^^^^ (All) 

Eqn. (jAlip is the master formula for the signal. Given a predicted 3-pt correlation function and detector setup we 
can calculate the expected signal to noise. 
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Taking the Fourier transform of the window function, W{t,t' ,t"), and breaking up our total integral into chunks 
of length AT > 1//^ we can perform the time integrals^" to obtain the signal per chunk 

(5i23)m = I dhdh'dn"fdff^df'f"^df" {h^{Lm^'{f\^')h^"if",^')) 

•' A,A' ,A" 

xF^{n)Ff {n')Ff {n'')e^'''''f^-'''+f'^' ■'''+^'''^''-'''^w{fj' J''). (A12) 

In general, the signal is reduced by the factors F^{n) and exp[27ri(/f2 • xi + f'U' ■ xa + f"n" ■ X3)]. F^{Cl) < 1 
because the detectors generally do not possess isotropic beams in general and detectors are often not optimally aligned. 
The relative location exponent, exp[27ri(/Jl • xi + f'fl' ■ X2 + f'tl" ■ X3)] < 1, as detectors that are not co-located can 
destructively interfere. This is because detectors in different locations are measuring different parts of the wave. In 
the calculation of the 2-pt correlator, the angular integral over these factors is collectively called the overlap reduction 
function [5lj, r(/), and encodes the effect of the detectors on the correlation. In the case of the 3-pt there is no 
easy distillation of this information into a simple factor (like r(/)) because in Eqn. (|A11[) . the correlation function 
itself involves a triangle constraint in momentum space (5(/J7 + f'tt' + f'fl") which makes separation of the angular 
information difhcult. 

As noted above, localization means that the noise component of the signal above is zero. T he noise is then the 
variance of the above estimator A^i23 = 5123 — (5123). By definition {N) = 0, and we take N = \/{Nf^. Defining as 
usual ([4i,[5|) 

nait) ^ J dfna{f)e-^^'f , with {na{f)nt{f')) ^ ^SabSif - f')N!{f ), (A13) 
where the indices a, b label the detector then the noise is 

{N?23) = J dfdrdr'^m{f)^N,{frm{f"rmfj'j")\^- (ai4) 

Given a 3-pt correlator, Eqns (jA12|l and (|Al4p are the ingredients we need to estimate detection sensitivity. 

Evaluating the expression for the signal in Eqn. (|A12[) involves messy calculations and requires knowledge of 
the detector set-up, locations and orientations. Given a signal and detector configuration, an optimal filter can be 
designed. In principle, each different 3-pt correlation function will yield a different optimal filter. For our purposes, 
this is an unnecessary level of detail. Instead, we will make several simplifying assumptions in order to estimate the 
signal to noise ratio: 

• We assume isotropic detectors, so F^{(1,) = 1 V {tl,A). In general, these factors can be small to vanishing if 
our detectors are misaligned. 

• We assume co- located detectors, so exp(/r2 • Xi -I- f'Cl' ■ X2 + f'Cl" • X3) = 1 via momentum conservation. 

• We use an idealized noise spectrum, where the noise N{f) is an inverse tophat with minimum N{f^) in the 
domain (/* - A/, /* + A/). We then choose simple tophat filters W = T^i(/)W^2(/)W^3(/), such that Wi{f) = 
W2{f) = WsC/) = ^(/* + A/ — f)9{f — + A/), where /* is the optimal noise frequency, i.e. the filter has 
support for the frequency range of 2A/ around this optimal frequency. Clearly this filter is suboptimal - the 
optimal filter depends on both the exact noise spectra, the configuration of the detectors and the actual signal 
itself. 

Using these assumptions, we can simplify the integration of Eqn. (|A12p . The delta function imposes a triangle 
condition on the momenta, defining a plane in the 3 dimensional momentum space. The signal is assumed to be 
isotropic, and so without any loss of generality, we can pick to point in the z direction. Focusing on the triangles 
defined on the planes orthogonal to this direction (see fig. [5]), the integral sums over all the possible shapes of a single 
triangle. We write 

dndn'dn"T{f, /, f")S{fn + fn' + f'n") ^ d cos e^^Sn^Tif, /, /")-^5(/ - Vf'^ + + 2/'/" cos e^), (A15) 



The results of the time integrations are really truncated delta functions. We ignore this here for the purposes of estimation. 
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z 



k' 




FIG. 2: The 3 momenta must form a triangle via conservation. By assuming isotropy of both signal and detector configuration 
(as we have done here), the integral Eqn. (|A12p can be reduced to an integral over all possible shapes of a triangle as shown 
in this figure. The measure is then completed by rotating around the z axis, and then around all possible directions of its axis, 
for a total of Stt^. 



where 6*^ is the angle between f2' and tl" and /', /") is defined by 

(/.(/, ti)h{r, n'Mf", n")) = T{f, /', f")s{fn + fn' + f'n"). (ai6) 

The factor of Stt^ is the measure obtained by revolving the triangle around the z axis {2tt), and then over the solid 
angle spanned by Ct (47r), as allowed by isotropy. The tophat filters then constrain the shapes of triangles which 
contribute to the signal, hence the range of 9/^. It is easy to see that this implies 

cosO^ - 2(/,-A/)2 ' ^'''^ ''"'^^ ~ 2(/, + A/)2 ■ ^^^^> 

Eqn. (|A12p becomes 



{Si23)m= I dcos0^ fdff''df'f"^df"T{fJJ")—Sif - + + 2/7" COS 0^). (A18) 



p 



The simple noise model means the 3-pt noise is given by 

{NI23) = Nl{U)Nl{f,)Nl{f,){Mf. (A19) 

Before we press on to calculate the prospects for detection, compare the estimator Eqn. (jA18|) to the more familiar 
2-pt SNR. The signal to noise for the 3-pt scales as (A/)'^/^ while for the 2-pt at first approximation do not scale 
with (A/) - hence the 3-pt SNR is more susceptible to sampling rate effects. As long as we keep the sampling rate 
X :S> f^,, this effect is negligible. 



APPENDIX B: SGW FROM SCALAR SOURCES 

In this work, we follow [2]| , and present some basic results regarding the generation of a 3-pt correlation function 
of stochastic gravitational waves by cosmological scalar fields. The key results are Eqns (jBlSP and its linear field 
limit, Eqn. 

The results we present here are completely general, and apply to any process in which a scalar field evades the no-go 
theorem of Dufaux et. al. [l^. Examples of such processes are those considered above, preheating and global phase 
transitions. Since preheating is a highly non-linear, non-perturbative phase, numerical simulations and evaluation 
of Eqn. (IBISP are required for accurate results. This is outside the scope of this paper. During the early part of 
reheating and during global phase transitions, however, the fields remain linear and thus it is possible to calculate 



13 



analytically. In this limit we obtain Eqn. (jB17[) . We evaluate this expression approximately for the case of the global 
phase transition of [H,[23]- A similar result can be obtained for preheating using the technology of [H, HH, however, 
one must numerically simulate occupation numbers of the preheating fields, and we leave it for future work. 

As noted above, gravitational wave production by scalar fields in a non-inflationary regime is substantially different 
from gravitational wave production during inflation. Quantum effects are negligible and purely classical effects, e.g. 
the relaxation of a disordered field or the turbulent motion of large masses, lead to the emission of gravitational waves. 
In this case, gravitational waves are sourced by the transverse traceless part of the anisotropic stress Ily 



where the energy momentum tensor for the inhomogencous scalar fields is the usual 
The equation of motion for the gravitational waves in an expanding background with the above source is 



/i"(k,r)+ fc 



a 
a 



hij{k,T) = 167rGa3n^'^(k) 



(Bl) 



(B2) 



(B3) 



where a prime, ', denotes a derivative with respect to conformal time and hij = ahij is the comoving metric amplitude. 
The traceless, transverse part of the anisotropic stress is found using the transverse traceless projector 



n//(k) = o,,,,,„(k)n,„,(k) = 



Pal}i)PjM - -Fy(k)P,™(k) 



n,™(k) 



(B4) 



where the projector is Pij (k) — dij — kikj . Now, as pointed out by the only relevant part of the energy momentum 
tensor is the product of the spatial derivatives, so 



a^n-(k) 



Piik)p,„,{k)--p,,ik)PiM 



(27r)3 J (27r)3 



pV"0a(p)0a(p')'5(k-p-p') 



(B5) 



Now, since the gravitational waves originate from relatively short intervals of time, in this work we neglect the 
expansion of space, dropping the a" /a factor 



v'le,J+k'vke^J = l&TiG a{T)T'^^ {k) 

where /iij(k, r) = Ufc(T)eij(k). We can construct the Green's functions for this Eqn. [H, [52 

1 



G(r;r') = - sin(fc(r - r')) 
k 



(B6) 



(B7) 



then, we have 



hij{T,k) = 



WttG 



dr'sin(fc(r-r'))a(T')r,f (k,r'). 



(B8) 
(B9) 

Now, supposing that the source is turned off at some time Tf, then the gravitational waves becoming freely propagating, 
hence the solution becomes 



KJT,k) = WttG dr'cos(fc(r-T'))a(r')T"(k,r') 



hij (k, ' 



a(r)— a(r) 
where Aij{k) and i?y (k) are obtained from matching the solution at t — Tf, 

WttG 



a{T) 



A,,{k) 
By(k) 



167rG 



dr' cos[k{Tf - T')]a(T')7;-^^(k, t'), 
dT'sin[fc(r/-r')]a(r')TTT(k,r'). 



(BIO) 

(Bll) 
(B12) 
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<^(P, r') 



AAAA/ 

/ ,^(p-k,r') 

FIG. 3: Vertex for generation of gravitational waves from linear scalars. Wiggly lines denote gravitational wave Green's 
functions while dashed lines denote scalar propagators (stochastic averages). 





FIG. 4: Unlike inflation, SGW sourced by scalars (or any active sources) are generated at 1-loop order. Each vertex is of 
order (hdi(i>dj(fy^'^ where TT indicates the traceless transverse part of the spatial partial derivatives on i.e. it is simply 
gravitational bremsstrahlung. The power spectrum generated by such a process will be a 1-loop diagram with two external 
GW legs. 

Now, today for frequencies / 3> we can Fourier transform, taking oq = 1 we find, in real space 

/i,.,(x,i) = / /^d/ /"drJe"2''*^*e2^*'^^/(Ay(k)sin(k-x)+By(k)cos(k-x)). (B13) 

J — OO J 

Working in the co- located approximation, we can choose x = 0. We can now calculate the 3-pt function 

(/.,,(k,/)/.,,(k',/')/^fc.(k",/")) (B14) 
(16^G)3 rf 



kk'k" 



■ ' dTdT'dT"a{T) sin[fcr]a(r') sin[fc'r>(r") sin[fc"T"] (TTT(k, t)TI^ (k' , r')rTT(k", r")) 



which leads us to consider the unequal time correlation function of three copies of the transverse traceless energy 
momentum tensor; 



(7^7(k, T)T7,T(k', r')rTT(k", r")) = 0,,,mn{k)0,k,op{k')Ok^,,r{k" 



d^p f d^q f d^ 



(27r)3 J (27r)3 J (27r)3 

a{p)Mk - p)M^)Mk' " q)0c(s)0c(k" - S) 



PmPnQoQpSqSr 

(B15) 



where we have dropped terms which vanish (by transversality of (k)). 



1. Linear Field Limit 



In the limit that the scalar fields are well described by Gaussian statistics, we can then evaluate the 6-pt function 
of the scalar fields using Wick's theorem and the scalar field propagator [3 

{MKr)Mk',T')) = Fab{k,T,T')6{k + -k'). (B16) 

In this limit, one can understand the process in terms of the graph in figlH Using the rules, 

1. Draw all diagrams, and label each vertex with a time. 

2. A vertex as in fig. [3] gets a factor of -~16TrGOij^jyin{k)pmPn, where ©(k) is the transverse traceless projector 
defined in Eqn. (|B4)) 
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FIG. 5: The gravitational wave power spectrum generated by a scalar source is given by the 1-loop diagram. 

3. An external graviton (gravitational wave) line gets a factor of the Green's function solution Gfc(r, r') to Eqn. 
(|B7[) . where r is the time at which the diagram is being evaluated and r' is the time associated with the vertex. 



4. An internal scalar line gets a "propagator," which is technically a stochastic average: (0q(p, r)(/if,(p, t') — 
Fab{p, T, r') 

5. Conserve external momenta with an overall delta, S(^^ hi) 

6. Integrate over internal (loop) momenta 

7. Integrate over times associated with each vertex from the initial time (when the interaction begins) to the final 
time when the source is turned off. 

We can evaluate the diagrams. The three point function is, after fourier transforming into frequency domain and 
using the co-located detector approximation 

(/i,,(/, n)hjk{f', Cl')hk^{f", Cl")) (B17) 

= -^0,,,^n{k)OJk,op{k')Ok^,gr{k")6{^k + k' + k")^^^ ^ ' dTidT2dT3a{n) sin[fcTi] a(T2 ) sin[fc'T2]a(T3) sin[A;"T3] 

X J ^P>"(p>''(p-fc)«(p-fc)''F,fc(|k'+p|,T3,T2)F„fc(p,Tl,T2)F,e(|k-p|,ri,T3) 

+ {p - kfip - fc) Wi^fccdk" + p|, r3, r2)i^ab(|k - p|, ri, r2)F,e(p, n, r3) 



Here the factor for 4 comes from the equivalent diagrams and k = /fJ etc. The power spectrum can be calculated in 
the analogous manner. The power spectrum is generated by the diagram in fig. [51 We obtain 

/i67rG\^ 

{hij(k,T)hij{'k',T)) = 2(5(k + k')Ofci,mn(k) ( — ^1 J dri J dT2a(Ti) sin(fcTi)a(T2) sin(fcr2) 

Xrf t| J ^pf^plp-^p-Fabip, Tl, T2)F,fc(|p - k|, Ti, T2). (B18) 

This result is consistent with that quoted in [Tsj . 
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